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We consider the motion of and Josephson current through a mechanically oscillating supercon- 
ducting island asymmetrically embedded in a Josephson junction. The electromechanical coupling is 
provided by distance dependent tunneling rates between the electrodes and the island. The system 
asymmetry, resulting from the geometrical configuration, leads, for weak coupling, to an equation 
of the mechanical motion that reduces to the well-known Duffing equation. At zero bias voltage 
the island motion is determined by the homogenous Duffing equation that opens up two separate 
regions of solutions depending on the superconducting phases. The island either moves under influ- 
ence of an anharmonic single well potential, or is governed by a double well potential that allows 
for off-center oscillations. Under applied bias voltage the island equation of motion turns into a 
modified Duffing equation, with time dependent coefficients, that demonstrate both quasi periodic 
and chaotic behavior. 
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I. INTRODUCTION 

Nano electromechanical system (NEMS) resonators 
may now be micro fabricated precise enough that the 
effects of tunneling electron coupling to the mechanical 
system are measurable^. Such dynamical interactions 
between the charge carriers and vibrational modes in a 
mesoscopic system have been observed in si ngle electron 
tunneling to suspended carbon nano tubeiP^. Similar 
effects of vibron-electron coupling possibly explain dif- 
ferential conductance dips and peaks in molecular elec- 
tronics devices^^ and differential conductance steps in 
STM based inelastic tun neling spectroscopy on local vi- 
bration modes on surfacea^HSI. 

Nano scale resonator setups are interesting as fast high 
sensitivity detection devices^ for mas ^IHSlD charge^, 
force^ and displacement!^ and as mechanical systems 
reach the quantum limit implications for quantum infor- 
mation technology may be tremendous®"^!. 

The field has further evolved to include and explore 
superconducting NEMS. One investigative direction has 
been to couple nano mechanical resonators to a supercon- 
ducting Cooper pair box^J or a superconducting quan- 
tum interference device (SQUIDpsl, in order to probe 
and control superconducting qubits, as well as detect dis- 
placements near the quantum limit. Another course has 
been to study Josephson currents coupled to mechanical, 
or molecular, oscillators situated within the tunneling 
junctiorPsHini The oscillator then acts to shuttle Cooper 
pairs at resonant levels. 

In this paper we introduce a double Josephson junction 
that is asymmetric with respect to the oscillatory motion 
of a superconducting island. The dynamics of the system 
is captured as the mechanical motion is coupled to the 
electron tunneling. Apart from reproducing the expected 
equation of motion terms found in Ref. [SUthe asymmetry 
adds a nonlinear cubic term. At zero bias voltage the is- 
land motion is consequently controlled by the relative su- 
perconducting phases through the Duffing equation. This 



equation, thoroughly studied in mathematics^^, has 
received a lot of attention in NEMS research since non- 
linear restoring forces act on small scale resonators, see 
Ref. US] for a review. The Duffing equation also shows up 
in driven macro sc ale re sonators that are geometrically 
similar to our setup™ 

At finite bias voltages the equation of motion is mod- 
ified to include harmonically time dependent coefficients 
to both the linear and cubic term as well as a harmonic 
driving force. No studies have been published on this 
Duffing equation variant to our knowledge. At Joseph- 
son frequencies above and below the eigen frequency of 
the oscillating island regular and stable quasi periodic 
motion is found, wheres more resonant frequencies yield 
chaotic solutions. Chaos is an inherent property of the 
driven Duffing equation^. 

The importance of nonlinearities and the Duffing equa- 
tion in a NEMS aspect comes from a number of sug- 
gested and investigated applications. Weak signal ampli- 
fication with low noise levels based on system sensitivity 
near bifurcation points is one active subjecP^HH] Other 
novel experiments utilize buckled nano resonator beams 
that oscillate within the confines of a double well poten- 
tial, typical to the Duffing equation, either to produce 
mechanical quantized qubit states in the resonator by 
cooling^ or to construct mechanical memory bits that 
work under room temperature by controlling transitions 
between the potential welLp^l 

This paper is outlined as follows. In Sec. [il] a detailed 
description of the mechanical system is given. Here we 
also derive the Josephson tunneling currents dependent 
on island position as well as the equation of motion for 
the island as it couples to the tunneling Cooper pairs. In 
Sec. |III| the solutions to the island equation of motion in 
absence of bias voltage are presented together with nu- 
merical results for the tunneling currents under different 
conditions including zero bias voltage with varying su- 
perconducting phases and finite bias voltage. In Sec. [TV| 
we summarize our findings. 
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II. THEORY 

A. Description of model 

The nanomcchanical system considered comprise three 
superconducting electrodes out of which two are fixed 
perpendicular to each other. At the intersection where 
fixed left (L) and right (R) electrodes point a third mov- 
able island is suspended by a cantilever. In absence of 
electromechanical coupling the island is allowed to vi- 
brate in the direction of the left lead with a restoring 
force proportional to the distance from equilibrium. The 
setup is inherently asymmetric with respect to the mo- 
tion of the island and an illustration of the system is 
shown in Fig. [I] 



The first three terms of H separately give the electronic 
structure of the leads in terms of BCS Hamiltonians, c.f. 
Eq. ([2]). The superconductors couple through tunnel- 
ing term Ht- Here, c K and c\. annihilate and create an 
electron in electrode x — L, I , R with momentum k and 
spin a. Electrons in the left, island, and right leads are 
denoted by momentum p, k, and q, respectively. A x is 
the superconducting pairing potential in lead x- 

The tunneling response to island vibrations is mod- 
eled by distance dependent tunneling matrix elements 
T p k and T q k. For small vibrations we use the linear ap- 
proximation 



P k 



(4) 



where a is a positive coupling constant and is the 
tunneling rate to the island at its equilibrium position. 
The matrix element for tunneling between the island and 
the right lead is given by 
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where R is the equilibrium distance from the island to 
the right lead. 

By assuming low temperatures, T ~ 0.01 — IK, we can 
work with vibrational energies of the uncoupled island in 
the range luq ~ 10~ 6 — 10~ 3 cV, which is small in relation 
to the typical electron energy of 1 eV. 



Josephson current modulated by the island 
oscillation 



FIG. 1: Schematic picture of the mechanically and electron- 
ically coupled system. The superconducting island (SC I) is 
free to move as indicated by the arrows, while the left (SC 
L) and right (SC R) superconductors are rigid, mi and ki 
denote the island mass and spring constant respectively. A 
bias voltage V can be applied. 

The system forms a double Josephson junction and we 
assume the island displacement u to be small compared 
to the distance between the superconducting parts. Our 
aim is to describe the tunneling current as the electronic 
process couples to the mechanical motion of the island 
and we begin by addressing the electron Hamiltonian of 
the system, 

H = H L + H R + H I + H T , (1) 

where, 

H l ,r,i = e « c Lx c ^ + [ A l,RJ c U c - k i + H - c - 

(2) 

H T = Y^ T pk C U Ck ° + r qk c q<r c k°- + H c - ( 3 ) 
pker qk(T 



We derive the tunneling current, defined by I x (t) = 
—e(N x (t)}, where N x is the number operator, at junction 
x(= L, R) to the island by following Ref. HO] and obtain, 

I x (t)=2eRet e -^ t+t '\[A(t), A{t')]) 

J-00 (6) 

+ e -^x(t-t')([A(t),A^(t')])dt', 

where lo x = n x —fii defines the voltage drop between lead 
X and the island (p x and \i\ are the chemical potentials 
of lead x an d the island, respectively). The operators 
M l ) = Y, K ka ^kcL W^O): for k, e L,R and the time- 
dependence is defined by 

c Ka {t) = e^c^e-^ 

where K x = H x — (i x N x and Kj = Hi — mNi. 

In Eq. ^ the junction current is divided into two 
terms which describe different tunneling mechanisms. 
The second term accounts for the single electron tunnel- 
ing and will not be addressed further in this text. Our 
focus is here devoted to the first term, which describes 
the Josephson tunneling current. 
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We make use of the Bogoliubov-Valatin transforma- 
tion c K(T = it K 7 KO - - Wnlng^ where T) = ±1 differs in 
sign for spin up or spin down electrons, whereas u K and 



v K are the coherence factors satisfying \u K 



1 



and u* K v K 



/(2E K ), where <j> x is the supercon- 
ducting phase in lead %. Through the transformation, we 
define the quasi-particle energies 



E K = y' (e K - n x ) 2 
We can, thus, write 

|A x ||Aj' 



|A X | 2 - 



([A(t),A(t')}) = 



k;k(7 



T« k (t)T Kk (t') 



(8) 



(9) 



where t = t — t' and 



'x — •fi~ ( fx i O ur assumptions 
of small vibrational energies justifies the approximation 
T p k(t') = Tpk(f) — TTpk(t), which leads to the Josephson 
current /^(i) from the left lead to the island 



I L (t) = J L [l - au] 2 sm(u j^ L t + cj) L ) 

+ T i\l — au]aucos (uij^t + 4>l)- 
Here, the amplitudes 

y^| r (o),2l A x ||A/| 



(10) 



J x (eV) 



Kk I 



2-B K i?k 



1 



1 



T x (eV) 



eV + E K + £ k 
(o) |2 l A x ||Ai| 



(11a) 



1 



[eV + E K + E k ) 2 (eV 



(lib) 

define the tunneling between the fixed electrode \ — L, R 
and the island in absence and presence of the coupling to 
the vibrational mode, respectively. 

The tunneling from the right lead to the island is given 
by Eq. ( 10 1 after replacing the tunneling matrix element 
Eq. Q with Eq. i.e. 

2 



Jfl(t) =Jr [}~a 



\/i? 2 +u 2 - R 



sin (wj, fl t + 4> R ) 



auu 



\/R 2 + u 2 -R) 

cos (uj, R t + (f> R ). 



(12) 



C. Island motion modulated by the electron 
coupling 

In addition to the cantilever spring force acting on the 
island, electromechanical coupling contributes with a dy- 
namical force. To model this more complicated equation 



of motion we construct a Hamiltonian that include en- 
ergy terms Hjl and Hj_ R originating from coupling in 
each junction, 



Hose = -ffosc + H.J-L + Hj, R , 



(13) 



where H^X = ^~ — I- P is the island momentum, 

while mi denotes its mass, and ki is the cantilever spring 
constant. 

Hj^l and Hj R are constructed out of the requirement, 



9 9Hj T 

2e^— = Ij, 



which is fulfilled if, 

Hj,l ~[l - au} 2 (1 - cos (u JlL t + 4> L )) 
Ze 



(14) 



2erri] 



[1 - au]apsm (u)j iL t + cj> R ) 



(15) 



and 



Hj. R = J ( I n 



y/R 2 + u 2 - R 



x (1 - cos{uj J:R t + (j) R )) 
2em/ 



+ — — (l-a \VR 2 + u 2 - R 
2emj \ L 

sm(u} J>R t + 4> R ). 



(16) 



Vr 2 + u 2 



With a complete Hamiltonian the full island motion, u, 
is obtained by solving the Hamilton equations of motion, 



dH„ 



dH n 



dp 



P = 



du 



(17) 



In doing so we arrive at the following differential equa- 
tion, 



mjii + (7l + 7 H 2 )ii + (fcj + 7fl,i)u = i^L, 



(18) 



where, 

1L = 
1R.2 

lR,l 
and 



sin (u JtJj t + 4>l), 



(1 + aR) 



VR 2 + u 2 



[R 2 



x sin {u3j, R t + 4>r), 
J R a ( 1 + aR 



< v vr 2 + u 2 
\ ^Ur R + 1 ' (t,s u ' /,/, ' / " /,,) 1 



-JiQ 



(1 + aw) 



T L Uj,L 

2J L 



+ 1 ) cos (w JjL t + £ ) - 1 
(19) 
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One may note that equation ( 18 ) lacks a driving force 



term, Fr, that will be present if the angle between the 
right lead and the island motion differs from 90°. 

The central island equation of motion contains both 
time and nontrivial position dependence in its coeffi- 
cients. A more transparent equation is found in the weak 
coupling and low bias voltage limit. Under such condi- 
tions a is small and T x -C J x , so terms proportional to 
either r x , aT x or a 2 are dropped to enlighten the terms 
of greatest physical relevance. 

We also bear in mind that u/R -C 1 and keep only the 
second order Taylor expansions, 



1 



Vr 2 



1 u z 
~R ~ 2R? 



(20) 



etc. 



The coefficients ( 19 1 approximate to, 



1L 

F L (t) 



0, 

0, 

eR 
J L a 



(cos (uj,Rt + <f> R ) - 1) I 1 
(cos {uj. L t + 4> L ) - 1) 



2R 2 J 



,r \ , (21) 




FIG. 2: Zero bias voltage phase diagrams of four characteristic 
solution regions to the equation of motion (231 as it depends 



on the coefficients A and B when F = 0. B is aiways negative 
and only sofutions from leftmost quadrants are physical. 



and in defining 
1 



A(t) = — [ki + k D (cos(ujj fit + 4>r) ■ 
mj 

2miR l 



(22) 



where, ko = (Jna)/(eR), act as a dynamical spring con- 
stant, we end up at the equation of motion, 



u + A(t)u - B(t)u 3 = F L {t)/mi. 



(23) 



This is a Duffing equation modified by time dependent 
coefficients. It is only analytically solvable for zero bias 
voltage, through series expansions^, or by Jacobi's ellip- 
tic functional 



III. RESULTS AND DISCUSSION 



The dynamics of equation ( 23 1 directly effects the over- 
all Josephson current through its solutions. Due to the 
equations nonlinear nature we approach these numeri- 
cally in the general case and analytically for zero bias 
voltage. 



A. System under zero bias voltage 

Even at zero bias voltage, A(t),B(t), F L (t) — > A, £?, F, 
the island equation of motion has a rich variety of 
solutions if the superconducting phase differences are 



nonzero. Depending on the coefficients sign in equation 
(231 this Josephson effect opens up two distinct regions 
of solutions, depicted in the phase diagrams of Fig. [2] 
While four regions of solutions are obtainable mathemat- 
ically, only two are physical since B is non-positive. 

The force term F is non-negative and acts to shift the 
island motion away from the left lead, consequently low- 
ering the DC tunneling rate. The following analytic so- 
lutions apply to the 4>l — case where F = 0. Note that 
the current II = at all times under such conditions 
within the approximations made above. 

Also note that and Tr are zero when no bias volt- 
age is applied which means that the second term in the 
current expressions vanish. 

In terms of energy the island is confined by the poten- 
tial V(u) = ^u 2 - f u 4 which for A > 0, B < is a 
single well. Under these conditions the equation of mo- 
tion ( 23 ) has one singular point of center type and all 
phase trajectories are closed. The equation is satisfied 
by the solution, 



u = uocn(f2i, fc), 



(24) 



where Uq is the amplitude, cn(x,y) is the Jacobi elliptic 
cosine, Q — \J A — Buq, and k — ~ B/2 ■ Uq/CI. See 
upper left quadrant of Fig. [2] 

If A < 0, B < the island motion is governed by a 
double well potential which give rise to three singular 
points of which two are centers, corresponding to the 
double well bottoms, while the third is of saddle type 
centered between the two wells. All phase trajectories 
are closed and as is clear from Fig. [2] solutions exist that 
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circumfere either one of the two singular points of center 
type as well as solutions that enclose all three singular 
points. For A < 0, B < solutions to equation ( 23 ) can 
be written, 

±U dn(u} 1 t, fci) for -y/^^y < < J 



A>0, B<0 



A<0, B<0 



M cn(w 2 i, k 2 ) for Uq > 



2\A\_ 



(25) 

where uj 1 = ^/-B/2u , uj 2 = yJ-Bu^ - \A\, k\ = w 2 /Vi, 
k% = u\ /w| and dn is a Jacobi elliptic function. 

The upper solutions above correspond to trajectories 
enclosing cither one of the two singular points of center 
type. A sign change on the initial condition u , within the 
limit, gives rise to oscillations of equal frequency whose 
origin is separated by a distance 2y[Aj/ (-B) in real 
space. A schematic picture of the two solutions are de- 
picted in Fig. [3j The different solutions will not change 
the tunneling current between the island and right lead 
but the DC component of the tunneling current between 
the left lead and the island is clearly affected. At </>l =0, 
or in other words = 0, this effect is absent since equa- 
tion (10 1 yields zero current. For a small phase shift 



4>l 7^ 0, on the other hand, both solutions exist together 
with a non zero tunneling current 1^. The magnitude 
of the DC tunneling current difference between solutions 
confined to the two separate potential wells is less than 
AI L , DC < J L 4ay/\A\/\-B)sm<l) L . 

Figure [4] illustrates how phase shifts 0^ distorts solu- 
tion trajectories as well as moves the singular points. On 
the negative side of the origin the singular point of center 
type moves in positive direction, as F grows larger, while 
the singular point of saddle type moves in negative direc- 




FIG. 3: Schematic image of the island trapped in one of two 
potential wells present in the A < 0, B < case, at zero bias 
voltage and small phase shifts (f>L. For the indicated limits 
on «o the island oscillates with its center either to the right 
(above) or to the left (below). 




FIG. 4: Phase portraits of the two possible solution regions 
as they depend on the size of the force term, F. The centre 
in the A > 0, B < case shifts toward the right while the 
phase trajectories become more elliptical with a major axis 
parallel to the velocity. In the A < 0, B < case the leftmost 
two centers vanish with growing F while the rightmost center 
slowly shift towards the right. 



tion. The two points eventually merge, leaving only the 
singular point of centre type on the positive side. This 
point, on the other hand, slowly shifts to more positive 
values as 4>l increases. 

For A > 0, B < the singular point of center type 
shifts towards positive values as 4>l increases at the same 
time as the phase trajectories distorts toward an elliptical 
shape with major axis along the u direction. 

One may also note that while the island oscillations 
pass the origin the frequency of the motion induced In 
is double that of 1^ due to the system geometry. As 
soon as oscillations are restricted to either the positive or 
negative side the tunneling currents have equal frequency. 

To analyze the Josephson tunneling under coupling to 
the mechanical motion of the island, at zero bias voltage, 
the Fourier transform is taken for a fine mesh of varying 
phase shifts <pR and </>l. For cf>L = the analytical solu- 
tions are used while a numerical solver is utilized when 
4>l 7^ 0. As far as values goes the results are to be taken 
qualitatively even though realistic input parameters are 
used. 

The current parameters and the coupling constant are 
set tcP 1 J L = J R = 0.1 mA and a = 0.01 A~\ As equi- 
librium distance between the leads R — 10 A is taken, 
while an island mass of mj — 1 fg is used. 

Figure [5] depicts the Fourier transform of the tunneling 
current Ir as a function of its frequency and the phase 
shift <pfi for a given initial value u = 0.1 nm and me- 
chanical spring constant kj = 0.01 N/m. Unless stated 
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0.5 1 1.5 2 2.5 3 



FIG. 5: Fourier transform of the Josephson tunneling current 
In as a function of its frequency and the phase shift (j>R. 
ljo — 6.6- 10~ 8 eV is the island eigenfrequency and 4>l — at 
all times. The figure is composed of two images divided as the 
color bars indicate. Both color bars are graded in \iA. The 
lower half correspond to single well solutions where A > and 
B < while the upper half stems from solutions confined to 
one of the the double wells in the A < and B < case. The 
inset is an enlargement of a small area around i/>r/' k = 1/2. 



otherwise these are the input values used in calculations 
throughout the remainder of this paper. As a result the 
energy associated with the eigenfrequency of the uncou- 
pled island is ujq = 6.6 • 10~ 8 eV. All figures of Fourier 
transformed currents are given without DC component. 
As <\)r increases from to ixj2 solutions to the island 
equation of motion are restricted to the A > 0, B < re- 
gion, while phase values between 4>r = n/2 and ir result 
in solutions within the A < 0, B < region having tra- 
jectories that enclose the positive singular point of center 
type. Tunneling currents are close to singly harmonic in 
the < 4>n < 7r/2 region and with increasing phase shift 
the frequency drops until it approaches zero as A vanishes 
when 4>ji = tt/2. The value (f>n = tt/2 has no significance 
in it self but stems from the parameter input kj = kjj. 

Above = tt/2 the island motion has a more com- 
plicated shape which is reflected in the larger number of 
harmonics needed for its description. The current ampli- 
tude, given in fj,A 7 is noticeably two orders of magnitude 
larger compared to the amplitude of the bottom half arc. 
This is expected for solutions to the island motion with 
trajectories enclosing one of the singular points of center 
type when the initial value lies close to the origin for the 
A < 0, B < case in comparison to the A > 0, B < 
case. The AC current lies superimposed on top of the 
DC current and while solution trajectories depicted in 
the bottom left quadrant of Fig. [2] moves away from the 
origin all the way around the singular point to the right, 
trajectories in the upper left quadrant never go beyond 



the initial value. 

The inset of Fig. [5] depicts a small area around </> r = 
7r/2 and indicates that the transition between the two 
arc structures is not direct. As A goes from positive to 
negative the island motion phase portrait in the lower left 
quadrant of Fig. [2] builds up from the origin. The two 
singular points of center type divide from the the single 
center and leaves a saddle point behind. This can be seen 
as a discontinues step in the lower arc of the inset and 
solution trajectories now enclose all three singular points. 
With bigger <j>R values A becomes more negative and the 
lying eight shape of the phase portrait grows. Eventually 
the separatrix curve reaches the initial value Uq where 
the island nears the origin infinitesimally slow and the 
frequency drops to zero. Shortly after this point is passed 
the oscillation amplitude drops to zero as the initial value 
u and the bottom of the double well potential coincide. 

The differences in amplitude between solution regions 
even out when input values are changed to Uq = 0.4 nm 
and kj = 0.017 N/m, in accord with the discussion above. 
For these values the transitions of the Fig. [5] inset hap- 
pens over a larger phase shift range as well. All por- 
trayed in Fig. [6] 

With non zero phase shifts, <\>l 0, the force term F in 
equation p3| becomes finite positive which changes the 
island vibration signature in the tunneling current. Fig- 
ure [7] illustrates this with four consecutive images where 
(j>L increases for each image from left to right. 

First off the current amplitude from solutions in the 
A > 0, B < region varies with growing cf>L and sec- 
ondly the distinctive features of each solution region be- 
come less pronounced. The changes can be understood 




1/8 ■ " 0.2 



I ' ' ' 1 — I 

0.4 0.8 1.2 1.6 2 

to/w. 

FIG. 6: Fourier transform of Ir as in Fig. [5] with initial 
condition and spring constant changed to uq = 0.4 nm and 
fcj = 0.017 N/m, which leads to uj = 8.6 ■ 10" 8 eV. Currents 
from the additional set of island motion solutions, enclosing 
all three singular points in the A < 0, B < region, clearly 
show up above the discontinues step in the lower arc. 
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FIG. 7: The Fourier transform of the Josephson tunneling current Ir as a function of its frequency uji and phase shift 4>R 
depicted with increasing phase shifts 4>l from left to right. The amplitude scale is logarithmic on the form log 10 (l + Ir* 10 s ). 



with figure [4] in mind. In the first two images of Fig. [7] 
the stable center in the A > 0, B < region moves in 
the positive direction toward the initial value uq and the 
amplitude diminishes. In the following two images the 
stable center has passed u and amplitude gets bigger. A 
nonzero complicates the picture further since it causes 
<f>ji to also shift the stable center towards more positive 
values. This is evident in the cfy^ = 0.4 image of Fig. [7] 
where the stable center passes the initial value uq just 
below 4>r = 7r/4. 

The tunneling current frequency never goes to zero in 
either image of Fig. [7]even though the double well poten- 
tial governs the island motion shortly above 4>r = 7r/2, 
and upwards, in the = 0.1 case. The positive force 
term F shifts the 4>R value at which separation of the sin- 
gular points occur. When this happens the center point 
of the single well potential has already passed the ini- 
tial value Mo and the separatrix curve is never crossed. 
In the remaining three images with higher values the 
double well potential never form and no separatrix curve 
appears. 

All images in Fig. [7]seem to indicate that the tunneling 
current frequencies perfectly match at the transition from 
region A > 0, B < to A < 0, B < in contrast to the 
cf>L = case. The discontinuity in Fig. [5] is very small 
on the other hand and numerical noise makes it hard to 
distinguish such fine features. 



B. System under finite bias voltage 

A bias voltage ujjx = —^>.j.r > over the junction 
setup significantly changes the character of the tunnel- 
ing current. Most noticeably the island motion is no 
longer strictly periodic, but rather quasi periodic or even 
chaotic. The current amplitude also becomes some three 
orders of magnitude greater than in the zero bias case. 
This is easily understood as the usual Josephson factor in 
expressions (10) and (12) varies between —1 and 1 while 



order of 10 3 . All figures presented below are obtained 
with 4>i — 4>h = 0. 

At low bias voltage, wj/w < 0.24 (wj < 1.6-10 -8 eV), 
the island motion follows a regular pattern where it is 
quasi periodic with a high frequency, low amplitude, os- 
cillation superimposed on a higher amplitude vibration 
of frequency equal to uij. This causes small ripples on, as 
well as distorts, the tunneling current dominated by the 
Josephson factor. The typical case situation is shown in 
the bottom four images of Fig. [9j 

The current induced coefficients in the island equation 
of motion are zero at 



t = 



2-k/uj 



0,1,2, 



(26) 



the factor associated with island motion changes in the 



where the island vibration is harmonic. The cantilever 
spring dictates the motion for a few periods until the 
nonlinear and driving force contributions rapidly grow in 
sync with the Josephson AC. In such time intervals of 
strong nonlinearity the phase portrait implies a stable 
center markedly shifted away from the left SC lead. 

Between 0.24 < ujj/loq < 6.1 (1.6 • 10~ 8 cV < wj < 
4.0 • 10 -7 eV) the harmonic driving force and island mo- 
tion frequencies are comparable, but instead of simply 
resonating, the island vibrates in a chaotic fashion. This 
behavior comes as no surprise for such a strongly non- 
linear driven system. Chaotic solutions are a well stud- 
ied property of the ordinary driven Duffing equation but 
here we only go as far as to compare Poincare maps 
taken at loj/ujq = 0.48 and ujj/ujq = 24 to conclude 
that the low voltage map is compliant with a charac- 
teristic chaotic map while the higher voltage map has a 
clearly quasi periodic structure. The middle four images 
of Fig. [9] illustrates the island motion and Josephson 
tunneling current in the chaotic region at loj/ujq — 0.48 
{loj = 3.2 • 10" 8 eV). 

Above uij/uiq > 6.1 (uj > 4.0 • 10~ 7 eV) regularity 
in the island motion reappears, as the top four images 
of Fig. [9j taken at wj/w = 24 {uj = 1.6 ■ 10~ 6 eV), 
indicate. In this region the Josephson frequency is higher 
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FIG. 8: Poincare maps of the island motion taken at the 
t = + n • 2tt/u)j n = 0, 1, 2, . . . intersection and bias 
voltages Uj/ua = 24 (left) and ojj/ljo = 0.48 (right). 

than that of the major oscillatory island motion and the 
tunneling current is subject to a slow modulation. 

In contrast to the low bias voltage case, where the 
momentary phase portraits vary adiabatically with re- 
spect to the island vibrations, high bias voltage causes 
rapid changes in the time dependent equation of motion 
coefficients. The comparatively slow island is subject 
to a quickly deforming single well whose bottom shifts 
from the origin to a finite positive value with period 
T = 2-k/ujj. With increasing bias voltages the fine wave 
pattern in the solution trajectories seen in the top left 
image of Fig. [9] diminishes. 

IV. SUMMARY 

We have studied electron tunneling coupled to a me- 
chanical oscillator in a double Josephson junction. The 
geometry of the setup is asymmetric with respect to the 
mechanical motion which introduces a nonlinear term in 
the oscillator equation of motion. 

The Josephson tunneling current over one junction is 
modeled as linearly dependent on the oscillator displace- 
ment directly while the gap width of the second junction 
changes as a hypothenuse to the displacement. An imme- 
diate consequence is that the mechanically induced cur- 
rent frequency over the first junction is half that of the 
second junction if the oscillator passes its equilibrium po- 
sition and that the frequencies are equal if the vibrations 
are restricted to either the positive or negative side. 

In the uncoupled system the oscillator is taken to vi- 
brate harmonically with a linear restoring force. Cou- 



pling adds a linear, cubic and force term to the equation 
of motion - all time dependent at finite bias voltage. In 
the zero bias voltage limit the differential equation re- 
duces to the Duffing equation if the super conducting 
phases differ. 

The homogenous Duffing equation has two sets of phys- 
ical solution regions in our setup. One with a single well 
potential if the linear term is positive and one with a dou- 
ble well potential if the linear term is negative. Which of 
the double wells the oscillator is vibrating in is indistin- 
guishable by looking at the alternating tunneling current. 
The DC contributions will however vary between the two. 

Superconducting phase shifts associated with the junc- 
tion in line with the direction of oscillations govern the 
force term and can be manipulated to shift the potential 
well bottoms away from the rigid SC lead. The phase 
shift associated with the second junction controls the lin- 
ear and cubic term and a sweep through to tt reveal 
that single well solutions are obtained at low phase shifts 
while larger values produce a double well potential if the 
mechanical spring constant is chosen properly. 

At nonzero bias voltage we find three domains of so- 
lutions to the island equation of motion with their own 
characteristics. For very low bias voltage, such that the 
Josephson frequency is low compared to the vibration fre- 
quency, the island motion is quasi periodic which distorts 
the tunneling current by superimposing small ripples on 
the current. In the intermediate voltage span, where the 
Josephson and island frequencies are of the same order, 
the system turns chaotic and the tunneling current gets 
irregular distortions. For larger voltages, such that the 
Josephson frequency is much larger than the vibration 
frequency, the island motion is again quasi-periodic. The 
current is roughly harmonic, however, with a slow modu- 
lation superimposed arising from the mechanical motion 
of the island. 

The present study is based on theoretical assumption 
that should be within the realms of the state-of-the-art 
experimental capabilities. Detection of chaotic dynam- 
ics in nanoscale systems would be interesting from many 
perspectives and it is with great confidence we anticipate 
experimental verification of our proposal. 
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